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Bosons in the form of ultra cold alkali atoms can be confined to a one dimensional (Id) domain 
by the use of harmonic traps. This motivates the study of the ground state occupations Xi of 
effective single particle states <f>i, in the theoretical Id impenetrable Bose gas. Both the system on 
a circle and the harmonically trapped system are considered. The Xi and (f>i are the eigenvalues 
and eigenfunctions respectively of the one body density matrix. We present a detailed numerical 
and analytic study of this problem. Our main results are the explicit scaled forms of the density 
matrices, from which it is deduced that for fixed i the occupations Xi are asymptotically proportional 
to i/N in both the circular and harmonically trapped cases. 
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I. INTRODUCTION 



The Bose Gas, whose genesis we gratefully owe to Bose and Einstein, is an icon of statistical physics. The Bose- 
Einstein condensate (BEC) has found its very special place in nature as the fundamental mechanism underlying the 
properties of a very special system, superfluid helium. In the past decade, a new chapter has been opened with the 
experimental realization of ultracold atomic Bose gases |l], || || . 

As such, new and intriguing possibilities for the Bose gas have arisen. As these systems require harmonic traps, it 
is possible to fashion them to constrain the bosons in low dimensional configurations. In particular, Olshanii jj] has 
pointed out that a one dimensional (Id) configuration can be obtained and maintained by making an elongated trap 
where one of the harmonic trap frequencies, oj z , is much much less than the other frequency uj±. The experimental 
realization of a Id ultra cold atomic Bose gas has now made this a physical reality^, [q, 

This exciting prospect has now brought back into prominence the noted quantum mechanical ./V-body problem, the 
impenetrable Id Bose gas. This system was introduced by Girardeau Q and Lieb and Linige r |p| . It was solved by 




Lenard in his two classic works |Ttj, 11 and in the ensuing major works of Sutherland []12l O, |14| , Vadiya and Tracy 
n |l7| and Jimbo et al @. 

a lc? system of N bosons in box of length L interacting via a repulsive delta- function potential of strength, 
was studied. It was shown that in the limit ep/gn — ► the impenetrable Bose gas was obtained, where is the 
Fermi energy and n = N/L is the number density. This limit corresponds to the condition na\ — * where g = ft 2 jma\ 
and ai is the lei s-wave scattering lengthQ. To constrain the bosons to a Id configuration the condition ep /?uu±_ — > 
must be obtained. This limit corresponds to the condition nig — > where Iq = y (fi/mwi), the characteristic range 
of the harmonic trap. These and concomitant matters are discussed in [l^, ^|. Appendix |X] gives the ground state 
energy and cf values for the finite N and large N limit systems. 

The signature of the Id impenetrable Bose system is the set of particle occupation numbers, which correspond in 
the untrapped system to the momentum distribution of the bosons but which is quite different for the harmonically 
trapped system. 

Section [n] presents our comprehensive stu dy of the occupation numbers for the finite N and large N limit untrapped 
impenetrable boson systems. Section III with Appendix [b] presents our comprehensive study of the occupation 
numbers for the finite N and large ./V limit harmonically trapped boson systems. Section [V presents concluding 
remarks. 
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II. IMPENETRABLE BOSONS ON A CIRCLE 
A. Properties of the density matrix 

1. The wave function 

We are interested in the ground state of a one dimensional many body Bose system which satisfies the free-particle 
Schrodinger equation, and also a condition of impenetrability which prevents two bosons occupying the same point 
in space. This system of impenetrable bosons is confined to lie in a box of length L, and satisfy periodic boundary 
conditions, i.e. 

il)%{xx,...,x k + L,...x N ) = ip N (xi,...,x k ,:-x N ), 

k = l,2,...,N (1) 

which may be interpreted as confining the particles to move on a circle of circumference length L. 

The impenetrability condition states that the wave function for our N body system must vanish whenever two 
coordinates coincide, 

ip%(x 1 ,x 2 , ...,x N ) = for Xj=x k {j^k). (2) 

For spin-less point particles however this condition is equivalent to the Pauli exclusion principle. Hence for any fixed 
ordering of the particles, e.g. 

x\ < x 2 < ••• < xn, (3) 

there is no distinction between impenetrable bosons and free fermions^. Therefore for the ordering (|^), the ground 
state of the system of impenetrable bosons can be constructed from a Slater determinant of distinct single particle 
plane wave states, which has zero total momentum and minimum total energy, 

f (i Ar 7V!)" 1 / 2 det[e 2mfc ^/ L ] fc j=i,..jv Nodd 

^(x 1: x 2 ,...,x N ) = I (iAr ^ !rl/2det[e2m(fe+1/2) r/Ay i, ^ i ;; ( 7 1,/2 Neven (4) 

k k=-N/2,...,N/2 

= {L N N\)-^ 2 Y[ 2i^7r(afc-a; J )/L]. (5) 

l<j<k<N 

Ignoring irrelevant constant phase factors, we note that in the region of configuration space which we are considering, 
defined by condition (||), the function (^|) is non-negative. This property of non-negativity distinguishes the ground 
state in Bose systems. Imposing the condition that the wave function for a system of bosons be symmetric under 
particle interchange, Xj <-> Xk (j ^ k), we then deduce from that for general ordering of the particles, the 
properly normalized TV body wave function for our system is 

^(x u x 2 ,...,x N ) = (N\L N )-^ 2 Y[ 2\sm[TT(x k -x j )/L}\ (6) 

l<j<k<N 
l<j<k<N 

This result was first obtained by Girardeau ||. 

2. Toeplitz determinant form for the density matrix 

The one body density matrix for a system of N > 2 particles is obtained by integrating out (N — 1) degrees 
of freedom from the product iPn(xi, x 2 , xn-i, x)iPn(xi, x 2 , xn-i, x'). Since the periodic boundary conditions 
imply translational invariance the one body density matrix will be a function of one variable only. Specifically, for 
the wave function (]?]) we have 

p%(x~y) = N I dxi... dx N -it/j N (xi,X2,:.,x N --i_,x)ip N (xi,X2,..;X N -i,y) (8) 
J Jo 
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l<j"<fc<W-l 



e «i _ e »»|2 



(9) 



where t = 7t(.t — y)/L, and to obtain (g) we have made a change of variables and utilized the periodicity of the 
integrand. The normalization of (Q) is chosen so that It is convenient in to introduce a dimensionless modification of 



g N (t) := Lp N (x 



(10) 



By utilizing the following general result it is straight forward to express @N(t) as the determinant of a Toeplitz 
matrix pl| , a result first noted by Lenard p0| , 



2-7T 2-ir 



j d6 x ... f<W N f[M) [] |e w '-e w *| a = det 

n (=1 l<n<m<JV 



■ 2tt 



,iey-fc) 



lo 



(11) 



j,k=l,2,...,N 



The result ( |11| ) can been obtained by identifying the double product as the squared modulus of the Vandermondc 
determinant, which implies that the left hand side of (Ell) is 



27T 27T 



i f Mi... f de N f[m) E E z( p MQ)f[ 

q o i=i PeSivQeSjv i=i 

N 27r 

= lut E E e(P)e(Q)Jl f dB f(B) 



= ie,[p(z)-Q(z)] 



»e[p(i)-Q(0] 



det 



2tt 

d0f{0)e ie< j- k ) 



i=i: 



j,*=l,...JV 



(12) 
(13) 
(14) 



where S'jy is the group of all permutations on TV objects, and e(P) is the signature of the permutation P. In obtaining 
(H) from @ we have used the general result that 



JV 



.TV 



^? E E e ( p ) e (O)n°i'(0.O(0 = det [%-,fe] i)fc =i„ 
' psSn QeSjy z=i 

Applying (|ll| ) to (||) and explicitly computing the one dimensional integrals we arrive at|i(i| 

g N (t) = det [7j-fc(<)]j )fc= i ) ... ) jv_i 



with 7n(t) given by 



7«(i) = - I 



2tt 



cos(u) - cos (t) \e lun du 

4cos(i) sin(n|i|) 



7„(t) = 2<J n , cos(t) - <5 n> i - £ n) _i - 

2 / sin([n + l]|t|) sin([n-l]|t|) 
7T V (n+1) (n-1) 



(15) 



(16) 



(17) 



(18) 



The succinct expressions ( |l^) an d ( |l8| ) are well suited for our numerical investigations of the momentum distribution, 
to be discussed in Section [IBS. 

Expanding out the determinant ( |l6| ) we find explicitly 



02 (t) = -[(7r-2i)cos(i) + 2sin(i)], 

7T 



(19) 



03 (*) 



1 

2^2 



f 15 + 2(tt - 2f) 2 + 4 (2 + tt - 2<) (tt - 2 (1 + £)) cos(2<) + cos(4i)) 
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+ -^(7r-2i)sin(2t), 



(20) 



QA(t) = 



1 



18tt 3 
1 

18tt 3 
1 

18tt 3 



271 + 36(tt - 2t) (tt - 2t) cos(t) + 9 -33 + 4(?r - 2t) (tt - 2t) cos(3i) 



[(tt - 2*) (25 cos(5i) + cos(7i) + 12 (tt - 2t) (3 sin(t) + 11 sin(3i)))] 



[8 (28cos(2i) + 5cos(4<) - 105) sin 3 (t) 



(21) 



For larger JV, ex panding out the determinant for gjy(i) quickly becomes unwieldy. We use the expressions (]l9|)-(pl|) 
in Section II B 1 to find the corresponding exact expressions for the momentum distribution. 



3. A differential equation for Quit) 

We present in this section a characterization of Qn(P) via a second order, second degree differential equation. This 
is a consequence of the identification |2^] of QN+i(t) with a certain well studied average over the eigenvalue probability 
density function of the circular unitary ensemble of N x TV random matrices. It has recently been demonstrated 
p3| that for a given N such an average is expressible in terms of a specific r-function occurring in the Hamiltonian 
formulation of the sixth Painleve equation, Pyi, and thus can be characterized by the Jimbo-Miwa-Okamoto cr-form 
of this equation. Hence we obtain as a corollary an explicit differential equation characterization for our QN(t). 
Specifically, as a consequence of the detailed discussion in p2| it can be seen that the function 

(e 2it -1) d 

Mt)= { 2i } ^log QN(t), (22) 
satisfies the second order nonlinear differential equation 

(1 - e-^fi-^it) + <j>%{t)f = 16 (l + <j> N {t) + 1(1- e-™) ^(t)) 
W 3 " U N {t) + ^ N {t)\ 4>' N {t) + (N-1)(N + 1) ($ N (t) + - (1 - e- 2lt ) ^(i)) ) • 

(23) 



The expression ( |22j ) can be inverted for gjv'(t) to yield 

g N (t)=Nexpl ds 6 N (s). (24) 



2% 
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If we substitute u = e into (23), and define 



a N (u) = 4> N {t) = u(u - l)-^-\ogg N (t)\ e 2H =u: (25) 

du 

then the resulting equation is precisely of the Jimbo-Miwa-Okamoto cr-form for Pyi- 

It is interesting to note p2]| that if in ( ^2| ) we replace QNif) with the one body density matrix for the free Fermi 
system, gff(t), then this function is also a solution of the differential equation (|23|), albeit with a different boundary 
condition. The boundary condition required to complete the characterization of £jv(i) or gfj(t) in terms of ( p2f ) and 
(^3|) is provided by knowing the first few terms of their small t expansion. For Qfj{t) this is trivial since we simply 
have 

- ^ - " (* - ^ + { ^f^<' + <*■>) ■ 

To obtain the small t expansion for fAr(i) we utilize the result, originally due to Lenard (ll], ^2), that p%(x) can be 
written in terms of an expansion over integrals of the n-body density matrix of the free Fermi system, which we denote 
here by p^{x, X2, X3, x n ) = L^^ n+1 ' Qj^(t,t2,ta, ...,t n ), with the obvious extension of our notation tj = irxj/L. 
Recall that g^(t, £2, £3, i n +i) has the following determinantal expression in terms of g^(t) 

r,F(++„+„ i .A-A*+( _ g N(t) [^(*i)]j=2,3,...,n+l 



IQnV' ~ tfcjjfc=2,3,...,n+l L^WW — l k )\j,k=2,3,...,n+l / 
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In the above, j and k refer to the row and column elements respectively, so that for example [^(tj)]j=2,3,...,n+i is an 
n dimensional row vector and [Qpf(tj — tk)]j fe=2 3,... n+i is an n x n matrix. 

The one-body density matrix for the impenetrable Bose system then has the form O, E2j, 



f n - n Jo Jo 



t 

F 



dt n +iQ N (t, t 2 , t 3 , .... t n+ i). (28) 

We remark [22fl t hat the right hand side of (|2S|) is exactly proportional to the first Fredholm minor for the kernel 
( p6| ) . Using fl26|)-(|28D , it is straightforward to deduce that qn (t) has the small t behavior 

^^-i^+^iW.. (29) 

which then serves as the boundary condition characterizing QN(t) in terms of (|22|) and (|2 



4- Asymptotic expansions of gjv(i), and the thermodynamic limit 



One important application of the differential equation (J23J) , is that it provides a far more expedient way of deriving 
the corrections to the small t expansion ( p9[ ) than the does the expression (28). By substituting a small t power series 
ansatz for 4>N{t) into the differential equation (|23j), we obtain equations which define all but one of the coefficients. In 
particular the resulting equation for the coefficient of t 3 vanishes identically, and to fix this parameter we require the 
boundary condition ( |29| ) obtained from (p8|). However the differential equation then provides an extremely efficient 
way of obtaining a large number of higher order terms in the small t expansion. We found it straight forward to obtain 
the first twenty terms of 4>N(t) in this manner, although it would require far too much space to exhibit all these here. 
Using the relationship ( [22] ) we obtain the corresponding small t expansion of qn (t) , and include here terms up to and 
including t 9 , 



QN(t) 
Qn(0) 



(TV 2 - l) „ TV (TV 2 - l) „ (3TV 4 - 107V 2 + 7) , 

1 - '-t 2 H i J -t 3 + ± J -t 4 

6 9tt 

N (lliV 4 — 40iV 2 + 29) ^ 5 (3N 



360 

,; -21iV 4 + 49iV 2 -31) 



1350tt 



15120 



N (N 2 - l) (1837V 4 - 12107V 2 + 2227) 7 



(TV 



7938007T 
1) (57V 6 - 557V 4 



2397V 2 - 381) TV 2 (TV 2 - 4) (TV 2 



1814400 



TV (TV 2 - 1) (253TV 6 - 3017TV 4 + 13867TV 2 - 22863) 
71442000tt 



24300tt 2 



t 9 + 0(t w ). 



(30) 



Comparing this expansion with the expansion for g^(t) given in (p6|), we see that while the terms with an odd 
power of t in (^0|) arise purely from the Bose nature of the system and are not present in the expansion of gjy(t), 
the even terms in ( |30| ) up to and including i 6 , are precisely just the corresponding free Fermi terms from ([26]). The 
coefficient of t s however is seen to be a composed of two pieces, the first is simply the corresponding free Fermi term 
which is rational, where as the second is irrational and dependent on the Bose nature of the system. Precisely the 
same behavior is observed in the well known small x expansion in the thermodynamic limit |lq| . 

Our differential equation (|2^) and small t expansion ( |30| ) for the finite system can be seen to reduce to the corre- 
sponding results in the thermodynamic limit as follows. We define 



ov(«) = lim 4>n{s/N) = lim ^(e 2 * 8 ^) 

N—^oo N—+oo 



(31) 



and denote the density matrix for the system in the thermodynamic limit by p(x). Taking TV — ► oo in ( p4| ) and 
(p3|), and retaining only leading order terms, we respectively recover both the expression for p(x) in terms of ay(s), 
and the correct tr-form of the fifth Painleve equation characterizing <Jy(s), as discussed in the celebrated work of 
Jimbo et al fl8|| . Further, recalling that the thermodynamic limit corresponds to TV — > oo with the Fermi momentum 
kp = irN/L held finite, we can replace Nt with kpx in (30) and then noting that only the highest powers of TV in 
each coefficient of ( J30| ) will survive in this limit, we recover the corresponding small x asymptotic expansion for p(x) 
due to Vaidya and Tracy O]. 
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In their tour de force jl7|, Vaidya and Tracy were also able to systematically construct the large x asymptotic 
expansion for p(x), the leading term of which is 

A-lfcFxI-V*. (32 ) 

The corresponding result for the finite system is a good deal more subtle, as we deal with competing behavior of 
large x and large N. Starting with the Toeplitz form (|T6|), a delicate analysis due to Lenardj24|], rigorously justified 
by Widom |25f| , gives the surprisingly simple result 

1/2 

Q N (t):=Lp%(x)~Npoo\Nsin(k F ^)\ = Poo VN\sm(t)\- 1/2 (33) 

with 

Poo = G4 (3/ 2 ) = ne i/2 2 -i/s A -e ra 92418 (34) 
v2 

where G(z) is the Barnes G function]^], A w 1.2824271 is Glaisher's constant and we have noted the identity [^6[ |27[ | 

G(3/2) = tt^V^ 1 / 24 ^ 3 / 2 m 1.06922. (35) 



It is now known that (16) is a special case of a class of Toeplitz determinants with singular generating functions 
for which the asymptotics can be computed (see for e.g. f28|| ). We note too that in the thermodynamic limit the 
expression ( |33"| ) reduces to (|32|). It is still an open question however, as to how to obtain the correction terms to 
( p3| ) and thus generalize the expansion of p[x) of Vaidya and Tracy to the finite system. These correction terms are 
important for investigating the large N behavior of the momentum distribution. 



B. Momentum distribution 

1. Preliminaries 



Due to our choice of periodic boundary conditions the Fourier coefficients of p^(x), 

f L (— 2mrxi\ 1 f v 
c n (N) = dx p%(x)exp I I = — / dt g N (t) cos(2nt), (36) 



L ) Wo 

have the physical interpretation of being the expectation of the number of particles in the single particle state of 
momentum 2imh/L. In terms of these momentum state occupations, we have the Fourier expansion for QNif) given 
by 

oo 

QN(t)= c n (N)cos(2nt). (37) 

71— — OO 

Recalling our normalization convention, ^tv(O) = TV, a direct consequence of ( j37[) is that 

oo 

N = £ c n (N), (38) 

71 — — OO 
OO 

q n (tt/2) = ]T i-lT^N). (39) 



71 — — OO 



The expression ( J38| ) simply states that the sum over the occupations of all momentum modes for our system is equal 
to the total number of particles TV. It is of interest to compare the asymptotic occupation of the even modes verses 
that of the odd modes, for large N. By subtracting (|3^) from (p8|), and observing that ( p3| ) implies qn{-k/2) ~ \ ~N, 
we obtain 

c 2n+ i(N) = — —. (40) 
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Hence for large N the total occupation of the odd modes is half the total number of particles, and so in this limit the 
particles are evenly distributed among both the odd and even modes. 

It is straight forward to obtain an exact expression for c„(2) for all n by performing the integration ([56|) with 02 (i) 
given by (|l9|), to obtain 

16 

C ^ = (4n» -!?*» • (41) 

We pause here to mention a mathematically interesting observation regarding the structure of the Fourier series for 
Q%(t) at the special values t = 0, and t = tt/2. By definition we know qn(0) = N, and as discussed in p2l the exact 
value of g»Ar(7r/2) can be expressed in closed form for all N. For N — 2 the particular expressions ( f4l| ) and (|3^) yield 

16 16 16 4 

ftj(0) = -fff?i(2) = 2 and g 2 (7r/2) = -^E 2 {1) = -j arctan(l) = - (42) 



where E\(s) and Ei(s) are the celebrated series studied by Euler|29, 3C[] 

00 00 (_-\\k 

^.(^EiiT^^'-n*) .>i, aw-g^j, (43) 

which have the property that the sum of E\(s) for even s > 2, and the sum of Sa(s) for odd s > 1 are both rational 
multiples of ir s . 

It is possible to obtain exact closed form expressions of c n (N) for other small values of N, however the complexity 
of the corresponding qn (t) obviously prohibits us from taking this procedure too far, and indeed already by N = 4 
the expression for c„(4) is quite a complicated object. Explicitly, 

^2 *o,n + r, |n| < 1 

Cn(3) = < 1 6n 2 + 2 , (44) 

^ H < 2 + n 2 (, 2 -i)V' " 1 

64 (80n 4 - 8n 2 + 45) 
(16n 4 -40n 2 + 9)V 

18874368 (l6?i 2 (26n 2 (88ti 2 (2ln 2 - 137) + 15613) - 179455) - 55125) 
(16n 4 - 296n 2 + 1225) 2 (16n 4 - 40n 2 + 9) 4 tt 4 



c„(4) 



2. Large N asymptotics of c n (N) 



Utilizing the large N asymptotic form for £jv(i) given by (p3[), we can obtain the asymptotic behavior of c n (N) for 
N > n, 



c„{N) ~ — VN [ dtcos(2nt)sm- 1/2 (t). 
n Jo 



This integral is expressible in terms of the beta function |3l}| , -B(a;, y), and so 

Poo 2 V27T cos(n7r) 



Cn(N) 



7T 5(n + 3/4, -n + 3/4) 



Vtt T(rt + 3/4) 
Thus the leading order large N behavior of Co (JV) is 



co(iV) 



T 2 (3/4) 



N « 1.54269V./V 



iV -> 00. 



(46) 

(47) 
(48) 

(49) 



A consequence of our lack of correction terms to the asymptotic result (33 ) for £>jv(i), is that only the leading order 
large N behavior of c n (N) is available to us via this argument. In Section [IB 3 we use numerical results to conjecture 
the form of the correction terms to ( 48 ) . 



8 





r (N\lr^(N\ pvalnatp_H liQinp 1 

L-7-1, 1 i V J / L.(J 1 1 V J vl V Ql U CI L LI Uolllti I \J\J } • 


r ( N \ 1 r»( N\ pvalnafpH nsinc*- l r ~\'~)\ 


1 


1/3 = 0.333333... 


0.332708... 


2 


5/21 = 0.238095... 


0.238061... 


3 


15/77 = 0.194805... 


0.194799... 


5 


221/1463 = 0.151059... 


0.151059... 



TABLE I: Comparison of the values for c n (N)/co(N) obtained via ( jHc| ) versus (p2|). 



Manipulating (|4^) we find that the ratio of gamma functions [Q reduces to a simple rational product and thus we 
obtain the following very tidy expression for the large TV behavior of c n (N) in terms of that of cq(N) as follows 



CnjN) 

co (AO 



n 

i=i 



(g ~ 3) 



(50) 



If n also is asymptotically large, so that N ^> n 1, we can simplify things further still by applying Stirling's 
formula to (E8I) which yields 



Cn (iV) 
CnjN) 

co (AO 



1 



64?r 



+ 0(n~ J ) 



r 2 (3/4) i 

V2tt 



64n' 



+ 0(n" 3 ) 



(51) 
(52) 



It is important to observe that there are no corrections linear in n _1 in the expressions ( pl| ) and (|52]), which is a 
consequence of the fact that the first correction terms in the Stirling expansions of T(n + 1/4) and T(n + 3/4) are 
identical. The practical consequence of this is that (^) provides an extremely good approximation of (^0|) even for 
small n, in particular even for n = 1. To demonstrate this remarkable agreement we list in Table | the values of 
c n (N)/c a (N) evaluated via both @ and @, for n = 1, 2, 3, 5. 

Finally, we compare our results for the finite system when N n > with the continuum momentum distribution 
constructed using the large x expansion ([32]) for p{x) of Vaidya and Tracy 0. Denoting the momentum distribution 
in the thermodynamic limit by c(fc), we find the following small |fc| behaviour 



c(k) 



—ikz 



p(x) dx ~ p 



\k F x\- 1/2 dx = \ - 



k 



-1/2 



(53) 



Substituting k = 2irn/L, consistent with our periodic boundary conditions, into (|5 
( |5l"| ) . This is consistent since even with n 3> 1 , if N 3> n this implies that k <C kp 



results in the leading term of 



3. Numerical calculation ofco(N), ci(N) and C2(N) 



In order to gain further insight into the occupations of the low lying momentum modes, we have performed the 
Fourier integral (|3^) over the Toeplitz determinant (H^) numerically, for various values of N, for the cases co(N), 
ci(iV), and c 2 {N). The results are shown in Figs. 0, |and |, where the dots represent the result of performing the 
numerical integration and the line represents a fit of this data to the ansatz 



c n (N) ~ aiVN + a 2 + a 3 N~ 1/2 + a^N^ 1 + 



(54) 



We also list in Table t he exa ct values of c (N) for TV = 2, 3, 4, 5, 6, 7. 

As discussed in Section [I B 2 , analytic arguments lead to \fN leading order behavior for the large N asymptotics 
of c n (N). The coefficients of y/ N for co(N), a(N), and C2(N) given in our numerical fits below match very well with 
those derived by such arguments. 



c o (A0 « 1.54273 Vn - 0.5725 4 
Cl (N) w 0.514345 Vn - 0.5739 
c 2 (A0 ~ 0.367622 ViV - 0.5775 - 



0.003677 



' N 
0.01128 



'N 
0.02948 



N 



(55) 
(56) 
(57) 
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10 20 30 40 50 60 70 80 
N 

FIG. 1: The zeroth mode of the momentum distribution co(iV) vs N. 




FIG. 2: The first mode of the momentum distribution ci(iV) vs N. 



The fits (|55|)-(|57]) were constructed using only the data for N > 30 so as to obtain asymptotic information for large 
N, however as can be seen from the plots they match up very well with all the data. Indeed, the x 2 value for each 
fit is of the order 10~ 10 . We note that the coefficient of the iV -1 / 2 term in each of (|55|)-(|57|) is considerably smaller 
than the coefficient of the leading term or the constant. Also we remark that the constant term is very similar in 
all three fits, and that a very similar constant again appears in the fits of the occupation numbers for the system of 



harmonically trapped impenetrable bosons to be discussed in Section [I B 3 . This suggests that the first correction 



term to the the occupation of the low lying states as a function of N may well be a universal constant common to 
both of these systems. 



4- Large n asymptotics of c n (N) 



[IB 2 



we found the asymptotic 



The occupation numbers, c„(iV), depend on the two parameters n and N. In Section 
behaviour of c n (N) for N ^> n. In this section, by contrast, we seek an expression for the large n behaviour of c n (N), 
with n 3> N. Observing the large n behaviour of the exact results for c„(2), c„(3) and c n (4), given respectively by 
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FIG. 3: The second mode of the momentum distribution C2(N) vs TV. 



N 


Exact co(N) 


Numerical co(N) 


2 


lb 

1 


1.62114... 


3 


1 | Jb 


2.10645... 


4 




2.51766... 


5 


1 | |[ WM 4 2U ,^4bU 


2.88069... 


6 


iy3bUV848()bS3()HubU4J.^§f8?^ r a^MTy^biJybyb | 4144 


3.20923... 




7 


7 1" 4nfifi3fi43328raino™i 31 75?nnnn ■** 1 anno 


3.51155... 



TABLE II: Values of c (iV) for TV = 2, 3, 4, 5, 6, 7 



(pd|), (HI) and fl45|), we are led to develop the expansion 

Cn(N) = 



„4 



02 



(58) 



To obtain the coefficients ai, 02, .. we construct a small t expansion of pnU) by applying the Mellin transform technique 
to its Fourier series representation (^) with Fourier coefficients given by (p8|), and then compare this result with (|3C|). 
Beginning with (37) we insert the Mellin integral representation of cos(2ni) which yields 



1 />C+200 

Co(N) + — / ds2 1 - s r(s)cos(7rs/2) 5 (s)^ s , 

^1 J C— 20C 



g N (t) 



where we have defined 



ff(s) = ^£r^2 = ^ a . c(s + 2j + 2) 



< c < 1 



(59) 



(60) 



and in the last step we have swapped orders of summation and identified the series representation of the Riemann 
zeta function, £( s )? which is finite everywhere on our contour of integration in (|59"|). 

Closing the contour of ( |59| ) to the left results in a small t expansion. The poles of the integrand in the left half plane 
are all simple, and lie at s = 0, —2, —3, —4, —5, .... We note that the term arising from the residue at s — combines 
with Co(N) to produce f?jv(0). Since T(s) cos(?ts/2) is analytic when s is odd, the poles at s = —(2k+ 1), k > 1, arise 
purely from the ak term in the zeta function expansion of g(s). Hence the corresponding t 2k+1 term in the small t 
expansion of ^jv(i) w iU contain ak, s o that we can obtain the values of all the simply by considering the odd terms. 
We note that as discussed in Section II A 4 , it is precisely the odd terms in the small t expansion of QN(t) which arise 
purely from the Bose nature of the system and do not appear in the corresponding free Fermi expansion. The poles 
at even values of s arise from T(s), since g(s) is analytic then, and so the residues at these values will contain the 
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entire g(s) series, (|60|). Calculating the required residues then, we obtain the following expansion 
^jv(O) TV 3N 1 3N 157V 2 45iV 

8tT 7 4fl(-8) a 47T q . . 

a 3 i 7 + — -t 8 a 4 t 9 . (61) 



315(iV+l) 315iV 28357V 

Comparing fl6l] ) with (^) we deduce 

iV 2 (7V 2 -l) 1 TV 2 (iV 2 - l) (-29+ 117V 2 ) 1 



Cn(^V) 



12tt 2 n 4 360tt 2 n 6 

JV 2 (N 2 - 1) (2227 - 121CW 2 + 1837V 4 ) 1 
20160tt 2 rfi 
N 2 (TV 2 - 1) (-22863 + 13867iV 2 - 3017iV 4 + 253JV 6 ) 1 
100800tt 2 n^' 



(62) 



Going out to higher orders is a straight forward matter, but the resulting expressions obtained for become 
increasingly cumbersome. 

Feeding the explicit values of from (p2|) into g(s) and comparing the even terms in (|6l]) with the corresponding 
even terms in ( |30| ) we obtain a consistency check on the validity of our ansatz. 

We note that for the values N = 2, 3, 4, the large n expansion ( |62| ) recovers the large n expansions one obtains 
from the exact exact results (pi), (|4|) and (fl5|). Further, we can also recover from ( |62| ) the large k behaviour of 
the continuum momentum distribution, c(fc), just as we recovered the small k behaviour of c(k) from ( |5l| ) in Section 
II B 2. If we apply the following asymptotic result p2| 

oo 

p(x) ~ p(0)J2bs(k F x) s , x^O (63) 

s=0 

poo i 00 fu \ 2 ( s +!) 

I p(x) cos(kx) dx ^(-l) s+1 (2s + l)!6 2s+1 ( J , k ^ oo (64) 

to the small x asymptotic expansion for p[x) of Vaidya and Tracy ||l7f 

p(x) _ (k F x) 2 \k F x\ 3 (k F x) 4 ll|fc F a;| 5 (k F xf l22\k F x\ 7 
p(0) ~~ 6 + 9vr + 120 1350tt 5040 + 105tt7! 

' h) (k F x) s ~ non lf J k F x\ Q + 0((k F xr) (65) 



24 300tt 2 9! / v ' 98 000tt27 2 ' 



4 (k F \\ 88 fk F \\ 244 /^\ 8 4048^ f k F \ 10 
c(fc) ^3^lTj + 45^lTj + 105^ I T + 1575^ I Tj " (66) 



we deduce that 



Substituting k — 2wn[L, which implies N/n — 2k F /k, into ( p2[ ) and taking N — > oo retaining only leading order 
terms then recovers (|6q). 

We can thus summarise the asymptotic behaviour of the occupation of the nth single particle momentum state for 
a system of impenetrable bosons on a circle as follows: for large N and fixed n, c n (N) diverges as y/N , whilst it 
tends to zero like n~ A when n 3> N. As will become evident in Section [II, this summary also describes the behaviour 
of a system of harmonically trapped impenetrable bosons on a line. 
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III. IMPENETRABLE BOSONS TRAPPED IN A HARMONIC WELL 
A. The analytic structure of the density matrix 

1. The wave function 

The wave function for a system of N impenetrable bosons on a line confined by a harmonic potential is defined by 
the Hamiltonian (in reduced units) 



A? 



N 



Q2 

3 = 1 3 3 = 1 



(67) 



and the impenetrability condition discussed in Section (If A 1). The well known single particle eigenstates of ( ply have 
the form 



j = 0, 1, 2, 



(68) 



where Hj(x) is the jth Hermitc polynomial. The ground state wave function for the harmonically trapped non 
interacting Fermi system is obtain ed by forming the Slater determinant of the functions ( |68| ) with j = 0, 1, N — 1. 
Arguing then as in Section (II A 1) we obtain the ground state wave function for the impenetrable bosons by simply 
taking the modulus of this Slater determinant. 

The Vandermonde determinant formula states that 



det\pj-i(xk)]j,k=i,...N = Yl ( Xk _ x i) 

l<j,k<N 



(69) 



for any set {pj(x)} with Pj(x) a monic polynomial of degree j, i.e. a polynomial of degree j for which the coefficient 
of is 1. In particular, ( p9| ) is true for {pj{x)} — {2~ j Hj(x)} and {Pj{x)} — {x j }. Applying ( |69| ) to the Slater 
determinant over ( |68| ) we see that the ground state wave function for N harmonically trapped impenetrable bosons 
can be expressed as 



1 N 

ipo(x 1 ,X2,-,x N ) = — w T\e~ x " /2 IT \xj-x k \, 



(70) 



N fc= 



where 



fc=i 



JV-l 



l<j<k<N 



(C§) 2 = N\~[[ 2- m ^m\. 



(71) 



m— 



2. The one body density matrix 
The one body density matrix is defined as in (||) with the domain of integration now R, 

/•OO f*OG 

p"(x,y) = N dxi... / dx N -i%l)"(x ll ...x N - 1 ,x)^"(xi,...x N -i 1 y)- 



(72) 



and is now genuinely a function of two variables. We are again able to find a closed form expression for pjy(x, y) in 
terms of a determinant, analogous to the result on the circle, this time the determinant being of Hankel type rather 
than Toeplitz type. To make this identification it is advantageous to note the following general result, in analogy to 

o 



1 N 

sfl 



AT! 



dxi g(xi) (det[fj-i(xk)]j,k=i,...,N) = det 



dtg(t)f j -i(t)f k - 1 (t) 



(73) 



j,k=l,...,N 
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which can be obtained by simply expanding the determinants and recalling ( |l5| ) so that left hand side becomes 

M E E e ( P ) e W)IT/ dx i 5(^)/QW-i(^)/FW-i(^) (74) 



det 



(75) 



j,k=l,...N 



To obtain a closed form for p^(x,y) we proceed as follows. Consider xjv-i, from (|7Fj) and factor out 

the pieces that depend upon x, then apply (By) to obtain 



AT 



-x 2 /2-y 2 /2 



N-l r00 



n 

i=i 



dxi e x < \x — xi\\y — x;|(det 



(76) 



and hence fl22fl 



2 (j+fe)/2 



2ViVr(j)r(A;) 



b h k(x,y) 



(77) 



j,k=l,...N—l 



where in the last step use has been made of (|73J). The elements of the determinant have the following explicit form 



bj,k{x,y) 
bj,k(x,y) 



dt e"' |ic-i||y-t|* j 



+/c-2 



1 < i, A < iV - 1 



df e"' (x-t){y-t)t 



j+k-2 



2sgn(y-x) dt e~ % (x - t)(y - t)t r 



-k-2 



bj,k(x,y) = fj,k{x,y) 

- 2 sgn(y - x) [xy^ j+k -2 (x,y)-(x + y)fi j+ k-i(x t y) + fx j+k (x, y)] 



where for clarity we have introduced 



fjA x >y) 



dt e~ t2 (x -t)(y- t)t j 



+k-2 



j+k-1 
2 



xy + r 

r(^) (x + y) 



j+fc+1 

2 



j+k even 
j+k odd 



Pm(x,y) 



dt e~ x t 



(sgn(y)) m+1 fm + 1 



(sgn(x)) 



m+l 



m + 1 



(m + 1) 



1*1 1 



m + 3 



x m+l £ -x 

(m + 1) 



r I -1 m+ 3 2 
1F1 1, — - — ,X 



(78) 

(79) 
(80) 

(81) 
(82) 
(83) 

(84) 



and where 7 and \F\ respectively denote the incomplete gamma function and confluent hypergeometric function. 

There are some properties of p^(x, y) worthy of note. Firstly, as is clear from the form of (|70|), (x\, X2, ■ xn)\ 2 
is precisely the probability density function for the distribution of eigenvalues of the Gaussian Unitary Ensemble 
(GUE) of random matrix theory, and hence p^(x,y) can be interpreted as a certain average over the GUE [p2| . 
Secondly it is obvious that since the ground states of the noninteracting Fermi and impenetrable Bose systems differ 
only in the presence of the absolute value in the wave function of the latter, we see by comparing (jT^) with ((79|) that 
the density matrix for a harmonically trapped noninteracting Fermi system is given by (|7^) with bj^{x, y) replaced by 
fj,k(x,y). From this observation it is immediately clear that the densities for the Fermi and Bose cases are identical 
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since on setting x = y in ( |80| ) we see that bj,k{x, x) — fj.k{x, x). The global large N limit of the density in this case 
results in the Wigner semi circle lawj|3|, Q 

J2N I x~2~ 

p*(x,x) = Z—^l-—, \x\<V2N, (85) 

while for \x\ > y2N we have p%(x, x) = to leading order. An oscillatory correction to this leading global asymptotic 
form has been given in [ |35"| . Finally we note from the second equality in (|8^) that since the \F\ function is entire, 
Pn(x,u) is analytic everywhere in the finite {x, y) plane except along the diagonal y — x, where it has discontinuities 
in its first derivatives. Such observations are important in choosing a numerical quadrature method. 

3. Occupation numbers 
Quite generally, for a many body quantum system the eigenvalue equation 

p N {x 1 y)4> j (y)dy = \j4>j(x), j = 0,1,2,... (86) 



defines the natural orbitals, (j)j(x), which have the physical interpretation of being effective single particle states, and 
the eigenvalues Xj which are interpreted as the occupation numbers for these natural orbitals. When the system lies 
on a circle the periodicity implies that the natural orbitals are simply plane waves, and so the eigenvalues are given by 
the Fourier coefficients of pn(x — y) and hence the momentum distribution and the set of natural orbital occupations 
coincide. In the general case, in which pN(x,y) is not translationally invariant, this is no longer true. Recent work 
|p6[ |37j] has focused attention on the computation and analysis of the momentum distribution for harmonically trapped 
impenetrable bosons. A more demanding task is the investigation of the N dependence of the occupation of the lowest 
natural orbital, Ao, the fundamental quantity of interest in discussing BEC-like coherence. We undertake a numerical 



analysis of A q in Section III B , and then in Section HID we go on to discuss how a new analytic result, obtained in 



Section [II C , yields its large N scaling. 

B. Numerical investigation of the eigenvalues of p^(x,y) 



The integral equation (|8^) was solved numerically for 2 < N < 30 to obtain the corresponding values of Ao (N) and 
Xi(N). In the present section we carefully analyse this numerical data in order to determine the N dependence of 
Ao(-ZV) and Ai(A^). Before getting under way however, a few comments are in order. 

1. Perfunctory remarks 

Since our ultimate goal is to infer from the sequence of computed A ( N) and Ai (N) values how these quantities 
grow with increasing N, it is not only the accuracy and number of the computed A values that are obtained that 
is important, but equally important is how we infer the large N behaviour from these computed values. Naive 
approaches such as fitting the data to a log/log plot in order to obtain the exponent in Ao ~ N a |}8| prove to be of 
no use. A striking demonstration of this fact occurs if we attempt to apply this procedure to our numerical co(N) 
data for the system on a circle. In this case we know a priori by analytic arguments that the exponent is precisely 
1/2. However, fitting the cq(N) data with 2 < N < 10 via a log/log plot predicts an exponent of a = 0.64, and even 
using 2 < N < 70 only forces the exponent to drop to a = 0.60. A similar analysis for the harmonically trapped 
system produces analogous results. Clearly such an approach will only be useful for very large iV, and more suitable 
methods of analysis are required. 

The method we shall use to analyse the N dependence of the computed Ao and Ai values is motivated by series 
analysis techniques. We assume an ansatz for Xq(N) of the form 

X (N)=aNP + b+^, (87) 

and similarly for Ai(iV). We then fit our numerical data to this ansatz, varying the parameters p and x and seeking 
to minimise the x 2 value for the fit. We have already noted that the fits of the numerical data for the system on the 



circle to this form of ansatz are in very good agreement with the known analytic results. We shall see in Section III D 
that again our numerical fits to this ansatz agree remarkably with analytic results that we obtain via independent 
means. 
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2. Numerical Results 

To numerically solve 

p oo 

p 1 N{ x ^y) < t ) ]{y)dy = Aj^(x), j = 0,1,2, . 



we used the so called quadrature method |39[ which replaces the integral equation with an approximate matrix 
equation. Specifically, we factored out e~ v / 2 in (|T?1 ) and applied a Gauss-Hermite quadrature (with abscissae £i, 
weights Wi and 1 < i < Z) to e y ^ 2 p^(x, y)4>{y). Setting x = £j in (|8|) for each i, we thus obtain Z linear equations, 
which then defines a matrix approximation to our integral operator. This matrix, S, can be so chosen that 



2 2 2 N ~ 1 

\m = e- x ^ 4 ' x ^^ 4 ^UjpUJ^ det 

V7T1 (7V J 



2 U+k)/2 



2v^Vr(j)r(fc) 



j,k=l,...N-l 

l,m=l,2,...Z, (89) 

where a suitable transformation has been used to force S to be symmetric, since such a transformation makes the 
numerical computations more efficient. We computed the required 6j,fc(£z,£ m ) using the representation in terms of 
incomplete gamma functions ( pO] ) with (p4|). 

Due to the complexity of S, each of its elements being given in terms of an (iV — 1) dimensional determinant, we 
find that as Z and TV are increased round off error becomes an important obstacle. For TV — 30 convergence was 
down to two significant figures and by TV = 35 we found that 64 bit precision became insufficient and to obtain higher 
values of TV a multiprecision routine would be necessary. As TV approached TV — 30 the fits of the data to (|S7j ) became 
rather sensitive to the precision of the Xj used. Due to this we chose to fit only the data 2 < TV < 27 to the ansatz 
( |S7| ) since for this set we could be sure that the precision was at least three significant figures (actually considerably 
more for the lower TV values). With the computed Xj analysed according to (|7|) however, it turns out that a good 
deal of large TV asymptotic information is already present even at these relatively low values of TV. This will become 



evident in Section HID 



Systematically varying the parameters p and x in (p7|), we observed a dramatic dip of three to four orders of 
magnitude in the x 2 of the fit when p hit the critical value of p = 0.5. For p = 0.45 and p = 0.55 the best x 2 values as 
x was varied were of the order 10 -3 , and these x 2 values became larger as p was moved away from p = 0.5 on either 
side. At p = 0.5 however, the best x 2 value was of the order 10~ 7 , which occurred for x in the neighbourhood of 2/3. 
This very steep dip in x 2 precisely in the neighbourhood of p = 0.5 indicates that the exponent of the leading order 
behaviour of Ao with TV is in fact 1/2, just as is the case on the circle. A similar analysis performed for Ai yielded 
the same sharp dip in x 2 a t P — 0.5, again dropping from around 10~ 3 on either side of p = 0.5 to around 10 -7 at 
p = 0.5. The coefficients of the correction terms for Ai turn out to be identical, to two decimal places, to those for Ao, 
however it now appears that the best value of a; is in the neighbourhood of x = 4/3. Specifically, fitting the computed 
Ao and Ai values with 2 < N < 27 we obtained 

A (iV) = 1.43^ -0.56+-^J (90) 



Ai(iV) =0.61^77-0.56+^^. (91) 

with x 2 ~ 10" 7 . Aj ;ain we remark that the similarity of the constant in these fits with that for the case of the uniform 
system, (|||) and (|56|), is quite suggestive. 

In FIGS. fUH we plot the computed numerical values for Xq(N) and Ai(iV) for 2 < N < 30 together with the fits 

©'©■ ™ . " X* 

Recently Papenbrock |37J discussed the computation of the density matrix using a form similar to (|77|). Information 

regarding the leading TV behaviour of Ao was inferred indirectly from a numerical investigation of the momentum 
distribution, by way of an interesting scaling argument, rather than by analysing the Ao(TV) directly as we have done. 
This argument predicts that Ao oc ^fN for large N, in agreement with our analytical study to be presented below. 

In the next Section we construct the large TV asymptotic form for p^(x,y). This result allows the vTV scaling of 
Xj to be put on a firm analytic footing as a straightforward corollary. 



10 




10 



15 
N 



20 



25 



30 



FIG. 4: Comparison of the computed values of Xo(N), with the fit of this data to the series ansatz, (f)o|). 
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FIG. 5: Comparison of the computed values of Ai(iV), with the fit of this data to the series ansatz, (f)o|). 



C. Global asymptotics from a log-gas picture 

We would like to derive the analogue of the asymptotic formula ( |33| ) for the impenetrable Bose gas in a harmonic 
trap. Now the latter is given as a determinant by (0). To our knowledge unlike the situation for the Toeplitz 
determinant (|l6|), there is no known general theorem giving the asymptotic form of determinants of the type J77|), or 
equivalently multidimensional integrals of the form (frgh. (If each factor \x — xi\\y — xi\ in (|7^ ) is replaced by e - ^*- 1 for 
g analytic, such general asymptotics are known pc[|). However, as we will now demonstrate, a Coulomb gas argument 
which can be used to give a heuristic derivation of ( j33|) does carry over to the case of the density matrix for the 
impenetrable Bose gas in a harmonic trap, thus allowing the global asymptotics of the latter to be predicted. 



Let us first revise the derivation of the asymptotic form (133), slightly modified relative to the presentation in 




so it is best adapted to application in the harmonic trap case. The starting point is the multidimensional integral 
formula (||) for QN{t). Put 
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/■27T /-27T n 

x / <»!.../ dSnTTle** 1 — e^ 1 1 291 |e^ 2 - e ie '\ 2q2 
Jo Jo , =1 

x | e «*_e <fl *| a (92) 

l<j'<fc<n 

This is the classical configuration integral for n mobile unit charges, and two fixed impurity charges (at angles </>i, 
4>2 and of charge qi, qz) on the unit circle interacting via the logarithmic potential — q a qb log \e l8a — e l6b \ at inverse 
temperature (3 = 2. In terms of this notation, and that of (pi), we can write 



r+ s 2nN Zg_ 1 ((0,l/2),(2t,l/2)) 
^(*)=|l_ e »|i /a ^ (93) 

where 

Z^:=Z^((.,0),(-,0)) = (2^) Ar 7V! (94) 

Now for N large and t fixed the impurity charges are effectively separated by a macroscopic amount (in units of the 
mean interparticle separation) so we expect the factorization 

Z^_ 1 ((0,l/2),(2t,l/2)) ^((0, 1/2)) Zg_ 1 ((2t, 1/2)) 

yC ~ yC yC 

N "N-l/2 "N-l/2 



' ^((0,1/2)) ^ 2 



(95) 



Now the reason for iV — 1/2 in the classical configuration integral of the denominators is to exactly balance the 
total charge in the numerator, as in the original expression (p3|). Of course to evaluate Z^_ 1 ^ 2 we must analytically 
continue off the integers but this is immediate from the exact evaluation in (p4|). Furthermore, we know from (^) that 

^((0,1/2) = ( ) [ G[N + 1 /2) ) 

~ {27r) N ^ 2 N~ 1 / 4 e~ N N N {G(3/2)) 2 (96) 
while it follows from (|94|) and Stirling's formula that 

Z^_ 1/2 = (2tt) n ~ 1 / 2 T{N + 1/2) - (2Tr) N N N e~ N (97) 



^((0,1/2)) (G(3/2)) 2 



and so 



Substituting ( |98| ) and (|9(j) in ( |95| ) then substituting the result in (^3|) reclaims ( p3[ ) 
This log-gas argument readily generalises to provide the global asymptotics for 



(98) 



g H (X, Y) := (2N) 1 / 2 p N (V2NX, V2NY) (99) 
with — 1 < X, Y < 1 fixed and N — ► oo. Let us write 

Z£ m ((Xi,<ft), (*;,,«,)) := | Xl „ X2 |2< Zl92e - m9l x 1 V2 e - m92 x 2 2 /2 

/2 <-2 « 

da*... / dx„ TTe"" 1 ^ 72 !^! - x/| 29l |X 2 - a;; | 292 
-2 J-2 

x [] |x fc -x,| 2 (100) 

l<j<fc<n 

where m = n + qi + q 2 . This is the classical configuration integral for n mobile unit charges and two fixed impu- 
rity charges (of charge 91,(72 at X±,X% respectively) interacting on the interval [—2,2] via the logarithmic potential 
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-<la<lb log \x a — X{,\ at inverse temperature (3 = 2. The charges are also subject to the one body potential (mg Q /4)x 2 . 



meet 

(We work in the scaled interval [-2,2] as this is the convention used in f42j| , the results of which we will be using 

Uj n\ 



subsequently.) Now it is well known from the random matrix interpretation of | 2 (recall ]22[) that the support for 



the density is the interval [— y2N, y2N]. Doing this then changing variables shows 

2N M ((2X,l/2),(2Y,l/2)) 



Q{X,Y) 



\2X-2Y\V 2 



7 H 
N.N 



where 



Z n ,m — Zn,m(( ' )0), (• , 0)). 

Now analogous to (^||), for N large and — 2 < X, Y < 2 fixed we expect the factorisation 

Z%_ hN ((2X, 1/2)(2F, 1/2)) Z§_ 1N _ 1/2 ((2X, 1/2)) Z%_ 1N _ 1/2 ((2Y, 1/2)) 



7 H 
N.N 



7 H 

Zj N-l/2,N-l/2 



7 H 

Zj N-l/2,N-l/2 



But in the notation of Brezin and Hikami (Eq. (46) of p2[|) 

^-i, J v-i/2(( 2 ^ 1 /2)) _ (N — 1/2)! 



7 H 

Zj N-l/2,N-l/2 



M 



Mi(2X) 



and it is shown in J42J (see also [|43j) that 

Afi(2X) ~ (2 7 r)- 1 / 2 (2iV(l - X 2 ) 1 / 2 ) 1 / 4 (G(3/2)) 2 

where we have used the fact pl[ 

K-l 



n 

1=0 



= (G(K+l)f 

(K + l)\ G(2K + 1) 



(101) 



(102) 



(103) 



(104) 



Substituting ( |i04[) in p03|) and noting that (N - 1/2)1/ N\ ~ iV" 1 / 2 , we then obtain from ( |102|) and ( |10l|) the sought 
asymptotic formula 



h (xy , Arl/2 (G(3/2)) 4 (l-* 2 ) 1 ^!-^) 1 /* 

? 1 ' J 7T |X-F|V2 



(105) 



D. Occupation numbers and natural orbitals as N — > oo 



Using ( |105| ) we can determine the large TV behaviour of the occupation numbers and natural orbitals. Taking the 
large N limit of the integral equation (|88|), utilising ( |105| ) and (|9^) we obtain 



(g(3/2)) 4 r 1 dx (i-x 2 ) 1/8 (i-n l/8 r m _ A, 

where we have defined the scaled limiting large N natural orbitals 

(2A0 1/4 <M*)-^-PO. 



(106) 



(107) 



This scaling for <f>Ax) arises since for large N the density matrix p^{x,y) has support \—y/2N, y/2 N], and so the 
integral equa tion ( p8| ) tells us that 4>j(x) must also, hence normalising 4>j{x) on this inte rval implies (107). 

From ( |106| ), the large N scaling of A., is immediate. Since both the left hand side of (|l06| ) and <fj(Y) are independent 
of N, we must have that for large N 

(108) 



Xj cx VN. 



We note that this argument has regarded j as a fixed parameter. For j » W we expect the behaviour of Xj to be 
similar to that of Cj(N) given by (|6l), since very excited states should not feel the presence of the trap. 
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FIG. 6: Comparison of the limiting large N natural orbital (po(X), denoted by the dashed line, with the the scaled natural 
orbital V2N(f>o(x) for N = 25. 



Hence from a knowledge of the asymptotic behaviour of p%{x, y) we have demonstrated not only that Ao ~ \^N, 
but that the other low lying modes must also display this N dependence. We can go one step further however and 
numerically solve ( |106| ) to obtain the prop ortionality constants missing from (108). 

Factoring out the exact constants, (106) becomes 



/ 



i 



IX-YW 2 



3 (1 _ y2)l/8' 



where 



We solved (109) using a similar method to that discussed in Section III B 2 , to obtain 

A = 3.438... Ai = 1.47... A 2 = 1.13... 
and hence, using [G(3/2)] 4 /7r = 0.4160... we find that in the large N limit 



A = 1.430VJV, 
Ai = 0.6lVN, 
A 2 = 0.47^. 



(109) 

(110) 
(111) 

(112) 
(113) 
(114) 



The results ( |112| ) and (113) are to be compared with the numerical fits ( |90| ) and ( |9l| ) for Ao and Ai. The agreement 
is remarkable. 

The limiting scaled natural orbitals are also of interest. FIGS. show <pj(X), together with (y/2N) 1/4 (f>j(x) 

computed from (§8|) with N = 25, for k = 0, 1, 2. It is clear that already by N = 25 the eigenf unctions of the large iV 
equation (106) provide a good approximation to the finite N results, only really differing outside the support. The 
cause of the sharp decrease as X — > ±1 can be understood by noting that (109) implies that <fj(X) vanishes like 
(1 - X 2 ) 1/S in this limit. 

We remark that ipj(X) is an even (odd) function when j is even (odd). For reference, we list the computed values 
of the central maxima and the area under the curve for <po(X) and ip2(X). Whilst the area under if2(X) is certainly 
smaller than that under ipo(X), it is not so much smaller as to be negligible. 



</?o(0) = 0.853... 



l¥ > 2 (0) = 1.20. 



J dX <p (X) = 1.38... J dX ip 2 {X) = 0.21.. 



(115) 



(116) 
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FIG. 7: Comparison of the limiting large N natural orbital ipi(X), denoted by the dashed line, with the the scaled natural 
orbital V2N(f>i(x) for N = 25. 



9- 




FIG. 8: Comparison of the limiting large TV natural orbital <p2(X), denoted by the dashed line, with the the scaled natural 
orbital y/Wfoix) for TV = 25. 



IV. CONCLUSION 



We have given an exhaustive study of the particle occupation numbers for the Id impenetrable boson systems. Our 
feature results are presented in Section III for the harmonically trapped systems. These results are of such strikingly 
similarity to those for the untrapped boson systems presented in Section || as to make all the results collectively 
invaluable. 

We commend the challenge to the experimentalists of measuring these signature quantities. 
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APPENDIX A: SUMMARY OF GROUND STATE ENERGIES OF SYSTEMS SUBJECT TO VARIOUS 

BOUNDARY CONDITIONS 



We list here for ease of reference the ground state energies per particle, Uq/N, and the Fermi energies, ep, for finite 
one dimensional systems of impenetrable bosons of mass m subject to periodic, Dirichlet, and Neumann boundary 
conditions and the analogous results in the thermodynamic limit, as well as for impenetrable bosons confined by a 
harmonic well. These results are identical to the analogous noninteracting Fermi results. 



System 


U /N 




Circle 


(N z - 1)tt^ /6mL z 


(N - l) z TY z h z /2mL z 


Dirichlet 


(N + 1)(N + l/2)n z h z /6mL z 


N z n z H z /2mL z 


Neumann 


(N - 1)(N + l/2)n z ti z /6mL z 


(N -l) z n z h z /2mL z 


Thermodynamic limit 


N z tt z ll z / QmL z 


N z TT z ri z /2mL z 


Harmonic well, V(z) — muj z z z /2 


Nhu z /2 


(N — l/2)huo z 



TABLE III: Summary of ground state energies per particle and Fermi energies for a systems of impenetrable bosons subject to 
periodic, Dirichlet and Neumann boundary conditions, and a harmonic well. 



APPENDIX B: TRIAL SERIES SOLUTION FOR A 

W e demonstrate here that it is possible to l earn a good deal about Ao by suitably approximating the integrand of 
(|109|). To begin with, simply taking Y = in (109) yields the following bound for Ao 



Ao = 
A < 2 



1 MX) (l-z 2 ) 1 / 8 

i MO) M 1/2 

1 (1_ 3,2)1/8 



dx 



11/2 



dx = 5(1/4,9/8) = 3.84108. 



(Bl) 



A very good approximation to the value of Ao can be obtained by inserting the formal Maclaurin expansion for 

MX), ' 



MX) = y>o(0)(l - AX 2 - BX 4 - ...) 



(B2) 



into ( |109|) and equating coefficients of powers of Y to obtain the unknown coefficients and Ao . Keeping only terms to 
order Y l and X 2 in the expansion of the integrand we thus obtain 



1+Y dt t fA l w 

[l-(A + -){t-Y) 2 ] + 
Vt o jq 



l-Y 



^[l-{A+\){t + Yf] 



\ [l-(A-l)Y 



which then yields, upon equating coefficients 



Ao 



Ao = 4-- 
5 





1 


3 




A - - 








8_ 


2 " 


h 2 





(B3) 

(B4) 
(B5) 



and the solution of these two simultaneous equations gives 

A = 0.57627... A = 3.43899... (B6) 
Carrying out this procedure again now keeping terms to fourth order we find that A and Ao shift only slightly 

A = 0.57682... 5 = 0.01914... (B7) 



A = 3.4378. 



(B8) 
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wh ich s hows that the p rocedure is quite stable. This value of Aq is in excellent agreement with our numerical solution 



of ( 109 ) given in (|lll[) . 
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